/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_PAM_sW_mC.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_PAM_sW_mC
{
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double n;
	double t;
	double tPrim;
	double U;
	double V;
	double h;
	double epsilonF;
	double beta;

	unsigned int L;
};

struct variables {
	double chiC;
	double chiCPrim;
	double chiCFup;
	double chiCFdown;
	double nF;
	double mu;
	double mC;
	double mF;
	double dF;

    double Energy;
 };

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.n = p.t = p.tPrim = p.U = p.V = p.h = p.epsilonF = p.beta = 0.; p.L = 0;
}

inline void clear(variables &v) {
	v.chiC = v.chiCPrim = v.chiCFup = v.chiCFdown = v.nF = v.mu = v.mC = v.mF = v.dF = 0;
	v.Energy = 0;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.chiC += q.chiC*d;
	v.chiCPrim += q.chiCPrim*d;
	v.chiCFup += q.chiCFup*d;
	v.chiCFdown += q.chiCFdown*d;
	v.nF += q.nF*d;
	v.mu += q.mu*d;
	v.mC += q.mC*d;
	v.mF += q.mF*d;
	v.dF += q.dF*d;
	v.Energy += q.Energy*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.chiC = q.chiC*d;
	v.chiCPrim = q.chiCPrim*d;
	v.chiCFup = q.chiCFup*d;
	v.chiCFdown = q.chiCFdown*d;
	v.nF = q.nF*d;
	v.mu = q.mu*d;
	v.mC = q.mC*d;
	v.mF = q.mF*d;
	v.dF = q.dF*d;
	v.Energy = q.Energy*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.n = q.n;
	v.t = q.t;
	v.tPrim = q.tPrim;
	v.U = q.U;
	v.V = q.V;
	v.h = q.h;
	v.epsilonF = q.epsilonF;
	v.beta = q.beta;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.t  p.tPrim  p.U  p.V  p.h  p.epsilonF  p.beta  p.L
	variables v; // variables  v.chiC  v.chiCPrim  v.chiCFup  v.chiCFdown  v.nF  v.mu  v.mC  v.mF  v.d  v.Energy
	double nF, dF, mF, mu; // contractions

	// Abbreviations - dependend on kx, ky
	double gk;
	double gkPrim;

public:
	double E_up[2];
	double dEdlchiC_up[2];
	double dEdlchiCPrim_up[2];
	double dEdlchiCFup_up[2];
	double dEdlchiCFdown_up[2];
	double dEdlnF_up[2];
	double dEdln_up[2];
	double dEdlmC_up[2];
	double dEdlmF_up[2];
	double dEdldF_up[2];

	double E_down[2];
	double dEdlchiC_down[2];
	double dEdlchiCPrim_down[2];
	double dEdlchiCFup_down[2];
	double dEdlchiCFdown_down[2];
	double dEdlnF_down[2];
	double dEdln_down[2];
	double dEdlmC_down[2];
	double dEdlmF_down[2];
	double dEdldF_down[2];

	double q_up, q_down;
	double dqdnF_up, dqdnF_down;
	double dqdmF_up, dqdmF_down;
	double dqddF_up, dqddF_down;

	double lchiC, lchiCPrim, lchiCFup, lchiCFdown, lnF, ln, lmC, lmF;

	Func(const parameters &p, const variables &v) {
		set(p, v);
		//refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fq(int s);
	double fdqdnF(int s);
	double fdqdmF(int s);
	double fdqddF(int s);

    // Energy
    double fE_up(double gk, double gkPrim, int s);
    double fdEdlchiC_up(double gk, double gkPrim, int s);
    double fdEdlchiCPrim_up(double gk, double gkPrim, int s);
    double fdEdlchiCFup_up(double gk, double gkPrim, int s);
    double fdEdlchiCFdown_up(double gk, double gkPrim, int s);
    double fdEdlnF_up(double gk, double gkPrim, int s);
    double fdEdln_up(double gk, double gkPrim, int s);
    double fdEdlmC_up(double gk, double gkPrim, int s);
    double fdEdlmF_up(double gk, double gkPrim, int s);

    double fE_down(double gk, double gkPrim, int s);
    double fdEdlchiC_down(double gk, double gkPrim, int s);
    double fdEdlchiCPrim_down(double gk, double gkPrim, int s);
    double fdEdlchiCFup_down(double gk, double gkPrim, int s);
    double fdEdlchiCFdown_down(double gk, double gkPrim, int s);
    double fdEdlnF_down(double gk, double gkPrim, int s);
    double fdEdln_down(double gk, double gkPrim, int s);
    double fdEdlmC_down(double gk, double gkPrim, int s);
    double fdEdlmF_down(double gk, double gkPrim, int s);


	// DEBUG
	void printLambda() {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (l)\n", lchiC, lchiCPrim, lchiCFup, lchiCFdown, ln, lnF, lmC, lmF);
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (q)\n", q_up, q_down, dqdnF_up, dqdnF_down, dqdmF_up, dqdmF_down, dqddF_up, dqddF_down);
	}

	// REFRESH
	void refresh() {
		nF = v.nF;
		dF = v.dF;
		mF = v.mF;
		mu = v.mu;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		q_up = fq(1); // 1 is up
		dqdnF_up = fdqdnF(1);
		dqdmF_up = fdqdmF(1);
		dqddF_up = fdqddF(1);
		q_down = fq(-1); // -1 is down
		dqdnF_down = fdqdnF(-1);
		dqdmF_down = fdqdmF(-1);
		dqddF_down = fdqddF(-1);

		if(DEEP_DEBUG) {
			printf("###-+- %10f %10f ||| ", q_up, fq(1));
		}

		// obliczanie poszczegolnych lambda:
		lchiC = -p.t;
		lchiCPrim = -p.tPrim;
		lchiCFup = -p.V*q_up;
		lchiCFdown = -p.V*q_down;
		lnF = -p.epsilonF - 2.*p.V*(v.chiCFup*dqdnF_up + v.chiCFdown*dqdnF_down);
		ln = 0;
		lmC = p.h;
		lmF = p.h - 2.*p.V*(v.chiCFup*dqdmF_up + v.chiCFdown*dqdmF_down);

		// zeby pozniej mozna bylo wypisac
		glob.qup =  q_up;
		glob.qdown = q_down;

		glob.lchiC = lchiC;
		glob.lchiCPrim = lchiCPrim;
		glob.lchiCFup = lchiCFup;
		glob.lchiCFdown = lchiCFdown;
		glob.lnF = lnF;
		glob.ln = ln;
		glob.lmC = lmC;
		glob.lmF = lmF;
	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		gk = 2.*(cos(kx) + cos(ky));
		gkPrim = 4.*cos(kx)*cos(ky);

		refresh();

		for(int i=0; i<2; i++) {
			E_up[i] = fE_up(gk, gkPrim, 1.-2.*i);
			dEdlchiC_up[i] = fdEdlchiC_up(gk, gkPrim, 1.-2.*i);
			dEdlchiCPrim_up[i] = fdEdlchiCPrim_up(gk, gkPrim, 1.-2.*i);
			dEdlchiCFup_up[i] = fdEdlchiCFup_up(gk, gkPrim, 1.-2.*i);
			dEdlchiCFdown_up[i] = fdEdlchiCFdown_up(gk, gkPrim, 1.-2.*i);
			dEdlnF_up[i] = fdEdlnF_up(gk, gkPrim, 1.-2.*i);
			dEdln_up[i] = fdEdln_up(gk, gkPrim, 1.-2.*i);
			dEdlmC_up[i] = fdEdlmC_up(gk, gkPrim, 1.-2.*i);
			dEdlmF_up[i] = fdEdlmF_up(gk, gkPrim, 1.-2.*i);

			E_down[i] = fE_down(gk, gkPrim, 1.-2.*i);
			dEdlchiC_down[i] = fdEdlchiC_down(gk, gkPrim, 1.-2.*i);
			dEdlchiCPrim_down[i] = fdEdlchiCPrim_down(gk, gkPrim, 1.-2.*i);
			dEdlchiCFup_down[i] = fdEdlchiCFup_down(gk, gkPrim, 1.-2.*i);
			dEdlchiCFdown_down[i] = fdEdlchiCFdown_down(gk, gkPrim, 1.-2.*i);
			dEdlnF_down[i] = fdEdlnF_down(gk, gkPrim, 1.-2.*i);
			dEdln_down[i] = fdEdln_down(gk, gkPrim, 1.-2.*i);
			dEdlmC_down[i] = fdEdlmC_down(gk, gkPrim, 1.-2.*i);
			dEdlmF_down[i] = fdEdlmF_down(gk, gkPrim, 1.-2.*i);
		}
	}
};


/////////////////
///////////////// FullSimplify
/////////////////

double Func::fq(int s) {
	return (1.414213562373095*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))))/Sqrt((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double Func::fdqdnF(int s) {
	return (1.414213562373095*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((0.25*dF)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + (0.25 + 0.75*Power(dF,2) - 0.5*nF - 0.25*mF*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		     0.7071067811865475*(1 - 1.*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))))/Power((nF + mF*s)*(1 - 0.5*(nF + mF*s)),1.5);
}

double Func::fdqdmF(int s) {
	return (1.414213562373095*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((-0.25*dF*s)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + (0.25*(1 + Power(dF,2) - nF)*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		     0.7071067811865475*s*(1 - 1.*nF - 1.*mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))))/Power((nF + mF*s)*(1 - 0.5*(nF + mF*s)),1.5);
}

double Func::fdqddF(int s) {
	return 0. + (1.414213562373095*(0. - Power(dF,2)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        (0. - 4*Power(dF,3) + dF*(-2. + 3.*nF + 1.*mF*s))/(2.*Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/Sqrt((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}
///////////////////////////////////////////////////////////////////////
double Func::fE_up(double gk, double gkPrim, int s) {
	return (-(gk*lchiC) - gkPrim*lchiCPrim - lmC - lmF - 2*ln - lnF - 2*mu - Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2))*s)/2.;
}

double Func::fE_down(double gk, double gkPrim, int s) {
	return (-(gk*lchiC) - gkPrim*lchiCPrim + lmC + lmF - 2*ln - lnF - 2*mu - Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2))*s)/2.;
}

double Func::fdEdlchiC_up(double gk, double gkPrim, int s) {
	return -(gk*(Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2)) + (gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF)*s))/
			   (2.*Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2)));
}

double Func::fdEdlchiC_down(double gk, double gkPrim, int s) {
	return -(gk*(Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2)) + (gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF)*s))/
			   (2.*Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2)));
}

double Func::fdEdlchiCPrim_up(double gk, double gkPrim, int s) {
	return -(gkPrim*(Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2)) + (gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF)*s))/
			   (2.*Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2)));
}

double Func::fdEdlchiCPrim_down(double gk, double gkPrim, int s) {
	return -(gkPrim*(Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2)) + (gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF)*s))/
			   (2.*Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2)));
}

double Func::fdEdlchiCFup_up(double gk, double gkPrim, int s) {
	return (-2*lchiCFup*s)/Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2));
}

double Func::fdEdlchiCFup_down(double gk, double gkPrim, int s) {
	return 0;
}

double Func::fdEdlchiCFdown_up(double gk, double gkPrim, int s) {
	return 0;
}

double Func::fdEdlchiCFdown_down(double gk, double gkPrim, int s) {
	return (-2*lchiCFdown*s)/Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2));
}

double Func::fdEdlnF_up(double gk, double gkPrim, int s) {
	return (-1 + ((gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF)*s)/Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2)))/2.;
}

double Func::fdEdlnF_down(double gk, double gkPrim, int s) {
	return (-1 + ((gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF)*s)/Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2)))/2.;
}

double Func::fdEdln_up(double gk, double gkPrim, int s) {
	return -1.;
}

double Func::fdEdln_down(double gk, double gkPrim, int s) {
	return -1.;
}

double Func::fdEdlmC_up(double gk, double gkPrim, int s) {
	return (-1 + ((-(gk*lchiC) - gkPrim*lchiCPrim - lmC + lmF + lnF)*s)/Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2)))/2.;
}

double Func::fdEdlmC_down(double gk, double gkPrim, int s) {
	return (1 + ((gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF)*s)/Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2)))/2.;
}

double Func::fdEdlmF_up(double gk, double gkPrim, int s) {
	return (-1 + ((gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF)*s)/Sqrt(4*Power(lchiCFup,2) + Power(gk*lchiC + gkPrim*lchiCPrim + lmC - lmF - lnF,2)))/2.;
}

double Func::fdEdlmF_down(double gk, double gkPrim, int s) {
	return (1 + ((-(gk*lchiC) - gkPrim*lchiCPrim + lmC - lmF + lnF)*s)/Sqrt(4*Power(lchiCFdown,2) + Power(gk*lchiC + gkPrim*lchiCPrim - lmC + lmF - lnF,2)))/2.;
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	//if(DEEP_DEBUG) { f->printLambda(); }

	double Etan_up[2] = { (1-tanh(p.beta*f->E_up[0]/2.))/2., (1-tanh(p.beta*f->E_up[1]/2.))/2.};
	double Etan_down[2] = { (1-tanh(p.beta*f->E_down[0]/2.))/2., (1-tanh(p.beta*f->E_down[1]/2.))/2.};

	for(int n=0; n<2; n++) {
		w.chiC += Etan_up[n]*f->dEdlchiC_up[n];
		w.chiCPrim += Etan_up[n]*f->dEdlchiCPrim_up[n];
		w.chiCFup += Etan_up[n]*f->dEdlchiCFup_up[n];
		w.chiCFdown += Etan_up[n]*f->dEdlchiCFdown_up[n];
		w.nF += Etan_up[n]*f->dEdlnF_up[n];
		w.mu += Etan_up[n];
		w.mC += Etan_up[n]*f->dEdlmC_up[n];
		w.mF += Etan_up[n]*f->dEdlmF_up[n];

		w.chiC += Etan_down[n]*f->dEdlchiC_down[n];
		w.chiCPrim += Etan_down[n]*f->dEdlchiCPrim_down[n];
		w.chiCFup += Etan_down[n]*f->dEdlchiCFup_down[n];
		w.chiCFdown += Etan_down[n]*f->dEdlchiCFdown_down[n];
		w.nF += Etan_down[n]*f->dEdlnF_down[n];
		w.mu += Etan_down[n];
		w.mC += Etan_down[n]*f->dEdlmC_down[n];
		w.mF += Etan_down[n]*f->dEdlmF_down[n];
		// w.d w tym miejscu jeszcze nie

        //if(f->E[n]>0) { w.Energy += -exp(-p.beta*f->E[n])/p.beta; } // basically 0
        //if(f->E[n]<0) { w.Energy += f->E[n]; }

        // it would be safer have following code - warning: 30 is chosen arbitrarily
		if(p.beta*f->E_up[n] > 30.) { w.Energy += -exp(-p.beta*f->E_up[n])/p.beta; } // basically 0
		else if(p.beta*f->E_up[n] < -30) { w.Energy += f->E_up[n]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E_up[n]))/p.beta; }

		if(p.beta*f->E_down[n] > 30.) { w.Energy += -exp(-p.beta*f->E_down[n])/p.beta; } // basically 0
		else if(p.beta*f->E_down[n] < -30) { w.Energy += f->E_down[n]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E_down[n]))/p.beta; }
	}


	if(DEEP_DEBUG) {
		//printf("###--- chiC=%.15e chiCPrim=%.15e chiCFup=%.15e chiCFdown=%.15e nF=%.15e mu=%.15e mC=%.15e mF=%.15e\n",
		//		w.chiC, w.chiCPrim, w.chiCFup, w.chiCFdown, w.nF, w.mu, w.mC, w.mF);
		printf("###--- %3d %3d %10f %10f \t | %10f ||", i, j, x, y, glob.qup);
		printf("###--- %.15e %.15e %.15e %.15e | %.15e %.15e\n", f->E_up[0], f->E_up[1], f->E_up[2], f->E_up[3], Etan_up[0], f->dEdlchiC_up[i]);
	}
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	unsigned int max = p.L/2; // attention: L should be even!

    // SMART WAY - ALSO OK
    // SERIAL CASE
   	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
   	// * edges (4 times)
   	for(unsigned int i=1; i<max; i++) {
   		onePoint(temp, i, i);
   		onePoint(temp, max, i);
   		onePoint(temp, i, 0);
   		//iii += 12;
   	}
      addTo(temp, w, 4);
      clear(temp);
      // * interior area (8 times) - SERIAL CASE
      for(unsigned int i=1; i<max; i++) {
   		for(unsigned int j=1; j<i; j++) {
   			onePoint(temp, i, j); //iii+=8;									//	sleep(1);
        	}
   	}
   	addTo(temp, w, 8);
   	clear(temp);
   	// * missing corners
   	onePoint(w, max, max);
   	onePoint(w, max, 0);
   	onePoint(w, 0, max);
   	onePoint(w, 0, 0);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.chiC /= L2;
	w.chiCPrim /= L2;
	w.chiCFup /= L2;
	w.chiCFdown /= L2;
	w.nF /= L2;
	w.mu /= L2;
	w.mC /= L2;
	w.mF /= L2;
    w.Energy /= L2;

	w.chiC += 8.*v.chiC;
	w.chiCPrim += 8.*v.chiCPrim;
	w.chiCFup += 2.*v.chiCFup;
	w.chiCFdown += 2.*v.chiCFdown;
	w.nF += v.nF;
	w.mu += -p.n;
	w.mC += v.mC;
	w.mF += v.mF;

	// calculate dd
	w.dF = 2.*v.dF*p.U + 2.*p.V*(v.chiCFup*f->dqddF_up + v.chiCFdown*f->dqddF_down);

	double W = 8.*( p.t*v.chiC + p.tPrim*v.chiCPrim) - p.h*(v.mC + v.mF) \
			 + p.epsilonF*v.nF + p.U*v.dF*v.dF + 2.*p.V*(v.chiCFup*f->q_up + v.chiCFdown*f->q_down);
	double Lambdy = 8.*(v.chiC*f->lchiC + v.chiCPrim*f->lchiCPrim) + 2.*(v.chiCFup*f->lchiCFup + v.chiCFdown*f->lchiCFdown) \
    		+ v.nF*f->lnF + p.n*f->ln + v.mC*f->lmC + v.mF*f->lmF;

    v.Energy = w.Energy + W + Lambdy;

/*
    printf("##%f | %f %f %f | %f ||", W, w.Energy, Lambdy, v.mu*p.n, v.Energy + v.mu*p.n);
    printf("##%f |\n", w.Energy+4.*v.mu*p.n);
 //   v.Energy = w.Energy \
 //               + 2.*p.V*(v.chiCFup*f->q_up + v.chiCFdown*f->q_down) \
 //               + 2.*p.V*(v.chiCFup*f->lchiCFup + v.chiCFdown*f->lchiCFdown) \
 //               + 8.*( p.t*v.chiC + p.tPrim*v.chiCPrim + v.chiC*f->lchiC + v.chiCPrim*f->lchiCPrim ) \
 //               + p.epsilonF*v.nF - p.h*(v.mC + v.mF) \
 //              + v.nF*f->lnF + p.n*f->ln + v.mC*f->lmC + v.mF*f->lmF + p.U*v.dF*v.dF;
*/


	if(DEEP_DEBUG) {
		printf("###+++ chiC=%.15e chiCPrim=%.15e chiCFup=%.15e chiCFdown=%.15e nF=%.15e mu=%.15e mC=%.15e mF=%.15e dF=%.15e\n",
				w.chiC, w.chiCPrim, w.chiCFup, w.chiCFdown, w.nF, w.mu, w.mC, w.mF, w.dF);
		//sleep(1);
	}

	//if(w.chiAB != w.chiAB) { w.chiAB = 1000; }
	//if(w.chiS != w.chiS) { w.chiS = 1000; }
	//if(w.chiT != w.chiT) { w.chiT = 1000; }
	//if(w.barDelta != w.barDelta) { w.barDelta = 1000; }
	//if(w.simDelta != w.simDelta) { w.simDelta = 1000; }
	//if(w.mu != w.mu) { w.mu = 1000; }
	//if(w.mAF != w.mAF) { w.mAF = 1000; }

    // results
    case_PAM_sW_mC::glob.Energy = v.Energy + v.mu*p.n;
	//glob.qt = f->fgt();
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //


int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.n = ((Par *) params)->n();
	p.t = ((Par *) params)->t();
	p.tPrim = ((Par *) params)->tPrim();
	p.U = ((Par *) params)->U();
	p.V = ((Par *) params)->V();
	p.h = ((Par *) params)->h();
	p.epsilonF = ((Par *) params)->epsilonF();
	p.beta = ((Par *) params)->beta();
	p.L = ((Par *) params)->L();

	variables v;
	v.chiC = gsl_vector_get (x, 0);
	v.chiCPrim = gsl_vector_get (x, 1);
	v.chiCFup = gsl_vector_get (x, 2);
	v.chiCFdown = gsl_vector_get (x, 3);
	v.nF = gsl_vector_get (x, 4);
	v.mu = gsl_vector_get (x, 5);
	v.mC = gsl_vector_get (x, 6);
	v.mF = gsl_vector_get (x, 7);
	v.dF = gsl_vector_get (x, 8);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7), gsl_vector_get (x, 8));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    E.run(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.chiC);
	gsl_vector_set (f, 1, w.chiCPrim);
	gsl_vector_set (f, 2, w.chiCFup);
	gsl_vector_set (f, 3, w.chiCFdown);
	gsl_vector_set (f, 4, w.nF);
	gsl_vector_set (f, 5, w.mu);
	gsl_vector_set (f, 6, w.mC);
	gsl_vector_set (f, 7, w.mF);
	gsl_vector_set (f, 8, w.dF);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7), gsl_vector_get (x, 8));
		printf("###  %24.15e# %25.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# | E = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			gsl_vector_get (f, 6), gsl_vector_get (f, 7), gsl_vector_get (f, 8), case_PAM_sW_mC::glob.Energy);
		printf("\n");
	}

    return GSL_SUCCESS;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: n=%f, t=%f, tPrim=%f, U=%f, V=%f, h=%f, epsilonF=%f, beta=%f, L=%d\n", n(), t(), tPrim(), U(), V(), h(), epsilonF(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "n   t   tPrim   U   V   h   epsilonF  beta   L\n");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "      chiC               chiCPrim                 chiCFup             chiCFdown                  nF                     mu              mC                    mF                   d");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_PAM_sW_mC::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.qup, glob.qdown, glob.lchiC, glob.lchiCPrim, glob.lchiCFup, glob.lchiCFdown, glob.lnF, glob.ln, glob.lmC, glob.lmF); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.qup, glob.qdown, glob.lchiC, glob.lchiCPrim, glob.lchiCFup, glob.lchiCFdown, glob.lnF, glob.ln, glob.lmC, glob.lmF); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                         E             SQRT(qup)            SQRT(qdown)             lchiC               lchiCPrim                 lchiCFup                 lchiCFdown             lnF                   ln                     lmC                  lmF\n");}

void example() {

	printf("\nEXAMPLE PAM_sW1\n");

    Par par(N_PAR);
    par.set_n(1.75);
    par.set_t(-1.);
    par.set_tPrim(0.25);
    par.set_U(12);
    par.set_V(-3);
    par.set_h(0.03);
    par.set_epsilonF(-2.);
    par.set_beta(500.);
    par.set_L(512);
    par.print_verboten();

    //-3.337653507840241
    //1.969706970284349e-01 2.167670420188718e-02 5.561015775038512e-02 6.897316869802131e-02 8.769845923441738e-01 6.919013389446085e-01 1.269263515538260e-02 1.414868211781073e-02 4.508986650056068e-02   -3.337653507840241     0.533274638256225     0.528492344170331     1.000000000000000    -0.250000000000000     1.599823914768675     1.585477032510994     0.771033272110589     0.000000000000000     0.030000000000000     0.019948130798299

    Var initV(N_VAR);
    initV.set_chiC(1.969706970284349e-01);
    initV.set_chiCPrim(2.167670420188718e-02);
    initV.set_chiCFup(5.561015775038512e-02);
    initV.set_chiCFdown(6.897316869802131e-02);
    initV.set_nF(8.769845923441738e-01);
    initV.set_mu(6.919013389446085e-01);
    initV.set_mC(1.269263515538260e-02);
    initV.set_mF(1.414868211781073e-02);
    initV.set_dF(4.508986650056068e-02);

    // Z 1D1t_v2: U=10, n=0.97
	//sdelta          barchiAB             barDelta                mu                     mAF                    d                      E
	//                 3.119786e-01         -3.807814e-23          4.999908e+00         -1.624410e-17          4.852439e-01             4.848583	// 0, silne d
	//                -3.918124e-01          2.251022e-01          3.970843e+00         -2.576902e-11          2.888477e-02            -0.368634	// D
	//                -3.859593e-01          1.206862e-01          3.804065e+00         -1.879365e-01          3.084875e-02            -0.363913	// D, m
	//                -4.051486e-01          1.318796e-15          4.143907e+00          7.110104e-13          3.659600e-02            -0.353778	// 0
    // E:
    //11.200000 0.999800    3.362494837648e-01    3.350551731900e-01    3.550356000366e+00    7.339835921625e-02    6.726979905711e-05   -0.230396672929960     0.001202721086514     3.830458471110253     0.001208435395048     1.957155709469804     0.162201431523197     0.160426865804345    -3.550517931441191     0.110931280330479


    printf("Starting point: "); initV.print(); printf("\n\n");
    Res res(N_VAR); res.print_label();
    RootGSL root(par, function, res);
    if(root.Solve(initV) == GSL_SUCCESS) { root.print(); }

    printf("Koniec\n");
}

void example_inf() {

    printf("#Start case_PAM_sW_mC - inf:\n");

    Par par(N_PAR);
    par.set_n(1.6);
    par.set_t(-1.);
    par.set_tPrim(0.);
    par.set_U(9);
    par.set_V(-2);
    par.set_h(0.0);
    par.set_epsilonF(-4.);
    par.set_beta(500.);
    par.set_L(256);
    par.print_verboten();

	// I'm choosing start point
	Var initV(N_VAR);

	// INFINITE
	Res res(N_VAR);
	RootGSL root(par, function, res);
	printf("#sdelta      "); res.print_label(); fflush(stdout);
	while(1) {
	    initV.set_chiC(frand()-.5);
	    initV.set_chiCPrim(frand()-.5);
	    initV.set_chiCFup(frand()-.5);
	    initV.set_chiCFdown(frand()-.5);
	    initV.set_nF(2.*frand());
	    initV.set_mu(8.*frand()-4.0);
	    initV.set_mC(2.*frand()-1.0);
	    initV.set_mF(2.*frand()-1.0);
	    initV.set_dF(frand());


		if(root.Solve(initV) == GSL_SUCCESS) {
			//printf("       "); root.print(); root.print(pFile); par.print_verboten(); fflush(pFile); fflush(stdout);
			printf("       "); root.print(); fflush(stdout);// exit(1);
		}
	}
}

} // end of namespace case2D2t_SGA_tPrim0
